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We use atomistic simulations to study the resonant acoustic modes and compare different calcu¬ 
lations of the acoustic mean-free path in amorphous systems with nanometric crystalline spherical 
inclusions. We show that the resonant acoustic properties are not a simple combination of the 
vibrations in the inclusions and in the amorphous matrix. The presence of the inclusion affects 
the transport properties mainly in the frequency range separating simple scattering from multiple 
scattering processes. However, propagation of acoustic wavepackets is spatially heterogeneous and 
shows that the amorphous/crystalline interface acts as a low energy pass filter slowing down the high 
kinetic energy motion whatever the vibration frequency. These heterogeneities cannot be catched 
by the mean free path, but still they must play an important role in thermal transport, thus raising 
the question of the correct modeling of thermal transport in composite systems. 

PACS numbers: 63.50.Lm, 63.20.Pw, 62.30.+d 


INTRODUCTION 


The understanding of the propagation of vibrational 
excitations in disordered materials represents still a chal¬ 
lenge for modern research. In dielectric systems these 
vibrations are responsible for the thermal transport, and 
thus the low value and the peculiar temperature depen¬ 
dence of the thermal conductivity in glasses [l| . Although 
many theories have been developed over the years, a com¬ 
plete understanding of the mechanisms responsible for 
the different behavior with respect to crystalline mate¬ 
rials has not been achieved. For longtime the so-called 
thermal anomalies (peak in the specific heat and plateau 
in the thermal conductivity around 10 K) have been as¬ 
cribed to the existence of glass-specific soft modes which 
pile up to give raise to a peak in the reduced density of 
states at energies in the 0.2-1 THz range: the Boson Peak 
(BP). Recent investigations suggest however that the BP 
does not arise from new modes but is the counterpart 
of the first Van Hove singularity of the corresponding 
crystalline system Q. Still, in the BP energy range vibra¬ 
tional modes appear to be strongly scattered, so that they 
are envelopes of different modes with a rapidly decreas¬ 
ing lifetime, which could be at the origin of the plateau 
in the temperature dependence of the thermal conduc¬ 
tivity. A classification between propagons (low energy 
plane waves modes), diffusons (reduced-life time modes 
beyond the BP energy) and locons (localized modes at 
high energy) has been proposed by Allen and Feldman 
and reported again quite recently by Parshin and Bel- 
tukov [3,0| • The understanding of thermal transport in 
glasses becomes an important issue also in view of the 
application of disordered materials in devices requiring 
very low thermal conductivities. Indeed, semi-conductive 
glasses have been recently investigated for thermoelec¬ 
tric applications [M, where the lowest possible ther¬ 


mal conductivity is required, associated with a good elec¬ 
tric conductivity. It is however difficult to get these two 
properties in a single system: for this reason composite 
materials made of alternating amorphous and crystalline 
systems can be promising. Indeed, super-lattices built 
from alternated crystalline and amorphous layers of sili¬ 
con have been recently proposed as a novel way to dete¬ 
riorate the thermal conductivity of silicon while keeping 
intact its electronic conductivity 12| . These applications 
raise the question of what happens to vibrational modes 
when they cross the interface between a disordered and 
an ordered material. Here we address this question com¬ 
paring the vibrational modes in amorphous silicon and 
in a composite material made of a crystalline inclusion 
embedded in an amorphous matrix. More specifically, 
we investigate the nature of the vibrational modes in the 
composite material and show that at high energy they dis¬ 
tinctly differ from the modes in the amorphous system. 
Interestingly, we find that the presence of the interface 
perturbs a travelling wave-packet slowing down the com¬ 
ponents with the highest kinetic energies, thus acting as 
a low-pass filter for thermal transport. 


In this article, we first focus on the vibrational reso¬ 
nant properties of the samples before investigating the 
effect of the crystalline inclusion on the propagation of 
mechanical waves. The article is organized as follows: 
the first section will present the numerical methods and 
samples. The second section will report the eigenmodes 
present in our samples through the computation of the 
vibrational density of states, participation ratio and vi¬ 
sualization of the eigenmodes shape. The last part will 
focus on the comparison between the mean free paths 
of mechanical waves obtained through different methods 
(dynamical structure factor versus wave propagation) in 
our silicon samples. 
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NUMERICAL MODEL AND TECHNICAL crystalline volume fraction, called CinASO. 

DETAILS 


In order to prepare the samples and compute the dy¬ 
namical propagation of waves at constant energy we per¬ 
formed classical molecular dynamics simulation. In par¬ 
allel, we have computed the resonant vibrations by exact 
diagonalization of the dynamical matrix of the system. 


Sample Preparation 


We studied an amorphous silicon system and systems 
composed of an amorphous matrix containing a spheri¬ 
cal crystalline inclusion. For this last configuration, sev¬ 
eral samples were generated. These systems consist of 
N Ki 11000 particles contained in cubic boxes of lengths 
Lx = Ly = Lz = 60 A. with periodic boundary condition 
(PBC). The interactions between particles are described 
by the Stillinger-Weber (SW) potential 0 frequently 
used to describe amorphous silicon. 

The method used to prepare the amorphous solid con¬ 
figuration is described in Ref 1^: a silicon crystal of 
262 144 particles was first melted at a temperature of 
3500 K, and cooled down to 10 K at constant volume 
at a quenching rate of 10^^ K/s. Then, the energy is 
minimized using a damped dynamics (see description be¬ 
low). A smaller cube, 60 A large containing Ri 11000 
particles is cut in this sample and a hole of the size of 
the wanted inclusion is carved in it. A pure crystalline 
sphere is generated through the repetition of a diamond¬ 
like primitive cell and is inserted in the spherical hole 
of the amorphous sample. At this point, the crystalline 
and amorphous parts are together, but the bounds at 
the amorphous/crystal boundary are not established, nor 
are the positions of the particles realistic with respect to 
the new system. The sample is thus annealed at 100 K 
during 100 ps, then, the potential energy of the system 
is minimized to insure the mechanical stability of the 
system. To realize this relaxation process the MD code 
is used to perform damped dynamics. Calculating the 
forces acting on each particles, we compute for each par¬ 
ticle the needed direction for minimizing its potential en¬ 
ergy. At each simulation step the particles are slowed 
down using damping forces opposed to their movements. 
This procedure can be written as a damped velocity Ver- 
let algorithm where the damping parameter, 7 , is selected 
to allow the particles to stop in a close potential energy 
basin. For our relaxation process we chose 7 = 0.61 THz. 
This relaxation process is performed until a potential en¬ 
ergy of 10 “® eV.A is reached, putting the system in a 
local minimum of potential energy with no particles close 
to plastic instabilities. 

Several samples were generated, changing the crys¬ 
talline inclusion size between 0 and 50% of the total vol¬ 
ume. Here we focus mainly on the sample with 30% 


Structural Characterization 

The investigated samples are listed in TableUwith their 
shear modulus, bulk modulus, density and sound speeds. 
We characterized their structure by computing the pair 
distribution function (PDF), as reported for the amor¬ 
phous and CinA30 sample in Fig. [T] For the amorphous 
sample, the distribution is characteristic of disordered 
materials with a local order shown by the sharp first 
neighbor peak followed by less ordered shells of neighbors. 
For both the amorphous and CinA30 sample, the nearest 
neighbor shell is maximum at 2.4 A and ends at 2.9 A in 
agreement with experimental results on amorphous sili¬ 
con 0. For CinA30, peaks at longer range can be seen, 
which can be related to the presence of the crystalline 
inclusion where periodicity imposes constant distance be¬ 
tween shells of neighbors. The first and second peaks are 
at the same position with and without inclusion. Indeed, 
even though the first sample is amorphous, the SW po¬ 
tential imposes the first neighbor distance and the very 
local organization. For both samples, the pair distribu¬ 
tion function shows a shoulder at 3 A. This shoulder 
appears when five fold coordinated atoms are present in 
the system 0 : this over-coordination seems to change 
the local order and shifts some particles toward this pref¬ 
erential position. It is interesting to notice that the peaks 
of the crystal in the CinA30 are at the same place that 
the ones of the pure crystal sample. In the inset of Fig. [T] 
we report the structure factor computed from the PDF. 
The maximum of the structure factor fixes the size of 
the pseudo Brillouin zone [l 6 | , which, for the amorphous 
sample, is given by the second peak at q = 3.8 A ^. 


RESONANT VIBRATIONAL PROPERTIES 


In order to understand the effects of the inclusion on 
the vibrational response we studied the eigenmodes of the 
systems. We started by computing the dynamical ma¬ 
trix (Hessian matrix) which, if anharmonic interactions 
between particles are neglected, can be defined as 0,0: 


^ct/3 _ 1 d^Epot 

dxfdx^ 


( 1 ) 


is thus the component of the dynamical matrix that 
links the displacement of the particle j in the direction 
/3 to the displacement of the particle i in the direction 
a through the potential energy variation induced by this 
displacement. If we assume the displacements in the sys¬ 
tem to be waves of the form 


uf = 


( 2 ) 
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Amorphous CinAlO CinA20 CinASO CinA40 CinASO 


Crist. Vol. fraction (%) 

0 

10 

20 

30 

40 

50 

Inclusion’s radius (A) 

0 

17.27 

21.77 

24.91 

27.42 

29.54 

Shear Modulus (GPa) 

29.0 

29.4 

29.5 

28.3 

29.9 

35.7 

Bulk Modulus {GPa) 

103.0 

105.0 

105.7 

105.8 

106.6 

108.1 

Density {g/cnP) 

2.31 

2.31 

2.31 

2.31 

2.31 

2.33 

Longitudinal sound speed {m/s) 

7831 

7901 

7924 

7883 

7963 

8175 

Transverse sound speed {m/s) 

3543 

3568 

3574 

3500 

3598 

3914 


TABLE I: Comparison of the structural properties of the amorphous sample and amorphous samples containing a crystalline 
volume fraction X from 0 to 50%, referred to as CinAX. 



FIG. 1: Pair distribution function and structure factor (inset) 
for the amorphous (red line), CinASO (blue line) samples and 
crystal (gray). The amplitude of the g(r) of the crystal has 
been decreased by a factor of 10 to indicate the position of 
the peaks in a pure crystalline sample. 

Newton’s law becomes 

~ (3) 

j,P 

y/rriirf and uP' are the eigenvectors and the eigenvalues 
of the dynamical matrix while w/27r are the eigenfrequen- 
cies of the system. Numerically, the dynamical matrix is 
computed in a discrete way, by taking the particles one 
by one, moving them 0.001 A from their equilibrium posi¬ 
tion and looking at the change of forces created on all the 
other particles, for all the directions. Then the exact di- 
agonalization of the dynamical matrix is done using the 
FEAST Eigenvalue Solver [l^, thus obtaining the iN 
eigenmodes of vibrations and eigenfrequencies of the sys¬ 
tem. For the computation of the quantities presented in 
the following (dynamical structure factor and mean free 
path), the interactions are described by the dynamical 
matrix, which is generated using the harmonic part of 
the SW potential. 



FIG. 2: Vibrational density of states for the crystal (gray 
line), amorphous (red line) and GinASO (blue line) samples. 
The purple line is the linear combination of the crystal and 
amorphous VDOS (See text for details). The inset represents 
VDOS!Lip and highlights the boson peak. 

Density of Vibrational Modes 

The vibrational density of states (VDOS) can be ob¬ 
tained by summing the number of eigenfrequencies ex¬ 
isting for each frequency interval 5uj = 0.11 THz. The 
result can be seen in Fig. [5] for the amorphous, CinASO 
and crystalline sample. The difference of the VDOS from 
the Debye’s model is represented as VDOS/uP in the in¬ 
set. An excess of frequency appears from 2 THz to 6 THz, 
corresponding to the so-called boson peak. 

The VDOS shows a similar shape for the amorphous 
and CinASO samples. The influence of the crystalline 
inclusion can be seen on the CinASO VDOS mainly at 
7.5 THz and from 11 THz to 17 THz. We plot in light 
blue the linear combination of the pure crystal and the 
pure amorphous VDOS meant to reproduce the CinASO 
VDOS, as done on nanocrystals by Stankov et al. [^ . 
This linear combination is defined as VDOSunear = 
e * VDOSc -I- (1 — e) * VDOSa- The values found for 
e in order for the VDOSunear to match at best the 
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FIG. 3: Weight e of the crystalline VDOS in a linear com¬ 
bination for matching the VDOS of systems containing from 
0% to 100% of crystalline inclusion. The line is the fit with a 
power law: f{x) = O.OOSSi^’^. 


FIG. 4: Vibrational density of states of the CinA30 (blue line) 
sample, and partial VDOS computed from the displacement 
of the atoms of the amorphous (red line) or crystalline (gray 
line) part of the sample. The pink line represents 70% of the 
VDOS of the amorphous sample. 


VDOScinAso are much smaller than the crystalline vo- 
lumic fraction (Fig. |31). It can also be noticed that even 
if the VDOSiinear and VDOScinAso match well, some 
discrepancies exist, in particular at frequencies where 
the crystal and the amorphous VDOS are very differ¬ 
ent. These two observations show that the VDOS of the 
modes of a sample with a crystalline inclusion is not a 
simple combination of the modes of the two independent 
crystalline and amorphous systems. To go further, we 
computed the partial VDOS for the amorphous and crys¬ 
talline part of the CinASO sample as: 

3N Nam 

gam{ui) = Y^ (4) 

2 = 1 j^am 


3N Ncr 

5 cr(w) = ^ ^ ||u;-f5(a;-a;*) (5) 

i— j^cr 

Fig. m shows that the partial VDOS of the amorphous 
part is quite similar to the VDOS of the amorphous sam¬ 
ple, reported here scaled for a better visibility, thus the 
crystalline inclusion does not have an influence on the 
modes of the amorphous part. On the other hand, the 
VDOS of the crystalline part is a very smoothed version 
of the one of the pure crystal. It is clear from this figure 
that the dips in the total VDOS are due to the crystalline 
part. This result raises the question of the spatial ex¬ 
tension of the eigenmodes and their possible localization 
at the interface. In order to investigate this point, we 
analyze more in details the eigenmodes in the following 
section. 


Vibrational Eigenmodes 


In order to quantify the spatial extent of the modes we 
compute the participation ratio (PR), which expresses 
the fraction of particles in the system taking part to the 
motion for a given mode. It is defined as: 


1 (Eti 


PR = 


2\2 


N 


2 ^ 2=1 


( 6 ) 


with Ui the displacement of the ith particle and N the 
number of particles in the system. The results obtained 
for the amorphous and CinASO samples are shown in 
Fig. m The general shape of the PR follows previously 
obtained results for amorphous silicon 0, 211. Among 
other differences, the PR of the CinASO sample is slightly 
shifted toward higher frequencies, reflecting the higher 
wave’s velocity found for this sample (See Table U). 

Five different regimes of vibration can be identified 
from the PR (See Fig. [S|); we report vibrational modes 
for these five frequency ranges in Fig. [5] and 0 for the 
amorphous and CinASO sample respectively. We show 
on these figures the displacement of the particles in a 
8 A thick layer of the sample containing the displace¬ 
ment of highest amplitude. The arrows direction and 
size represent the direction and amplitude of the parti¬ 
cles displacement. 

Regime (1): At zero frequency, PR=1 for the three 
translational modes corresponding to the global trans¬ 
lation of the system in the three directions of space. 
Up to 2.0 THz, there is a mix of high PR (delocalized) 
modes and low PR modes that can be identified as soft 
modes [l^. Plane waves are delocalized over the whole 
system leading to the high but wavelength dependent PR 
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FIG. 5: Participation ratio (PR) for the amorphous (red line) 
and CinASO (blue line) sample. The CinASO PR has been 
shifted along the y-axis for the sake of visibility. The dashed 
lines separate the five frequency regimes discussed in the text. 


seen at these frequencies. Soft modes are the superposi¬ 
tion of delocalized modes of long wavelength, and dis¬ 
placements involving only small selected groups of parti¬ 
cles. These modes are known to be the markers of immi¬ 
nent plastic deformation (their frequency tends to zero 
when the system approaches a plastic deformation) 
and it is known that their number is much larger in a 
silicon-like than in a Lennard-Jones system [ 2 ^. It is in¬ 
teresting to notice that the influence of the inclusion is 
not limited only to the interface. 

Regime (2): This range of frequency shows a broad 
area of high PR that reaches 55%. Very characteristic 
whirlpool-like structures can be observed here (Fig.[ 6 ]and 
Fig. Ill), which arise from the disordered arrangement of 
the particles in the amorphous matrix and have already 
been studied by one of us 22, 2^. For the CinASO sam¬ 
ple, the crystalline and the amorphous matrix move to¬ 
gether. The modes extend over the whole sample leading 
to a high PR. In this frequency range, the presence of 
an amorphous/crystal interface does not seem to prevent 
the propagation of the modes. This is probably related 
to the fact that both systems have large values of the 
VDOS in this frequency range. Recent works 2^ show 
that for this range of frequency the modes are almost 
entirely longitudinal. 

Regime (3): Here we observe a drop of participation 
ratio that reaches a local minimum of 35% for both the 
amorphous and CinASO samples, independent from the 
presence of the inclusion. The displacements are here 
localized but less spatially focused than soft modes. 

Regime (4): This region of the PR is characterized 
by two interesting drops around 12.5 and 15 THz for 
the CinASO sample. Precursors of these drops can be 
seen in the amorphous sample’s PR but they become 
more marked as the size of the inclusion increases. The 


largest displacements, that show where the modes are 
localized, are on groups of particles in the amorphous 
matrix. These modes are free to extend in the purely 
amorphous sample, where the PR is therefore high. How¬ 
ever, in the case of CinASO, the crystalline part stays very 
still, thus reducing the PR. It could be due to the fact 
that the crystal has fewer modes of vibration at these fre¬ 
quencies |4j the modes of vibration cannot extend in the 
crystalline part where they do not exist. It should be no¬ 
ticed that similar drops of PR have already been observed 
in amorphous silicon by Allen et al. and have been in¬ 
terpreted as ’’resonant modes” that are not localized but 
temporarily trapped in regions of peculiar coordination. 
For Allen et al. this interpretation is valid for the drops 
of PR at low and high frequency. We show here that 
the four different low PR areas can be explained by dif¬ 
ferent phenomena. More recently, similar drops of PR 
have been observed in nanocrystalline silicon where their 
origin was the confinement of these modes on grains in¬ 
terfaces or on grain interiors 2]J, as well as in Si/Ge core¬ 
shell nanowires where no interpretation was proposed 
and in silicon nanotubes where the low PR modes were 
localized on the surface of the material 26|. The origi¬ 
nal result here is that these modes are localized in the 
matrix rather than at the interface. Still, one could won¬ 
der what would be observed in a crystalline matrix with 
amorphous inclusion. 

Regime (5): For these frequencies, the PR shows a 
continuous decay to 0% up to 18 THz, corresponding 
to modes more and more localized as the frequency in¬ 
creases. As it can be seen in Figg. | 6 ] and [7] the localiza¬ 
tion is now organized in patches containing each several 
particles while the localized modes observed at lower fre¬ 
quencies were organized around small spots of localiza¬ 
tion. The high frequency localized modes can be seen as 
a manifestation of Anderson localization for phonons and 
have started to be studied in the framework of the multi¬ 
fractal theory applied to participation ratio The 

shape of these modes in the CAinSO sample is more diffi¬ 
cult to interpret because of the presence of optic modes 
in the crystalline part at these frequencies. 


MEAN FREE PATH 
Dynamical Structure Factor 

To further study the vibrational properties of our sys¬ 
tems, we computed the dynamical structure factor using 
the method described by Y. Beltukov et al. in Q. At 
the first time step of the MD run, once the structure 
is equilibrated, we impose to all atoms speeds randomly 
generated using a Maxwell-Boltzmann distribution cen¬ 
tered at 100.0 K with variance 0.2 K^. The forces be¬ 
tween particles are described by the dynamical matrix: 
only harmonic interactions are taken into account here, 





























6 



FIG. 6: Representation of different modes of vibration for the amorphous sample. The frequencies of the modes are from left 
to right 0.795 (Regime 1), 4.817 (Regime 2), 7.663 (Regime 3), 12.107 (Regime 4) and 17.451 THz (Regime 5). 



FIG. 7; Representation of different modes of vibration for the CinA30 sample. The frequencies of the modes are from left to 
right 0.816 (Regime 1), 4.903 (Regime 2), 7.803 (Regime 3), 12.205 (Regime 4) and 17.511 THz (Regime 5). 


and non-linear processes are neglected, thus the MD sim¬ 
ulation runs at constant energy (NVE) with a neglige- 
able temperature variation. The output of this run is a 
hie containing the positions of the particles at each time 
step. Using these data we can compute the dynamical 
structure factor, S{q,uj), for the system by performing a 
spatial and temporal Fourier transform on the displace¬ 
ment of all the particles for each timestep starting from 
the beginning of the MD run: 


S{q,uj) = 

2 ^ fT 

— ^ J u(rl,t) •m,exp(-igTl)exp(iwt)dt 


(7) 


where T is the total duration of the MD run and ffiq 
a polarization vector parallel or perpendicular to the 
wavevector q to have access separately to the longitu¬ 
dinal or transverse waves. Following this procedure, we 
obtained the longitudinal and transverse dispersion re¬ 
lations (Fig. [5]). A folding of the acoustic modes can 
be seen for both samples, centered at q = 1.9 A and 
q = 2.2 K for the longitudinal dispersion relation in 
the amorphous and CinASO sample respectively. These 
values are compatible with the pseudo-Brillouin zone bor¬ 
der, as determined by the most intense peak of the static 
structure factor, at Ri 3.8 A (inset of Fig. [T]). Generally 
speaking the dispersion relation of the CinASO sample is 
better defined than in the amorphous one. This is most 
striking for the transverse waves: for both samples there 
is a sudden broadening of the dispersion at ~ 3 THz, 
suggesting a transition from a propagative to a diffusive 
regime [^. However in CinASO newly propagative-like 
modes reappear between 10 and 11 THz, while above 


only non-dispersive modes are present, compatibly with 
the fact that only optic modes exist at these frequencies 
in the crystal. 

In order to get quantitative information from our re¬ 
sults, we have calculated the frequency and linewidths 
of acoustic modes using two procedures: one is to do 
a statistical analysis of the theoretical data; the second 
one is to artihcially create a pseudo-experimental spec¬ 
trum that we can thus analyze using standard fitting 
procedures. For both methods, the acoustic dispersion 
is obtained from the maximum of the dynamical struc¬ 
ture factor at each q (Fig. |9]). As for the width, different 
calculation methods can also be used on the theoretical 
data [i,[2i,[3^. 

For the first method, we compute the width with the 
standard deviation formula, as: 


Auj{q) 


r 

Jo 


'' (oj— < UJ >q)'^S{io) du) 


1/2 




( 8 ) 


For the second method, in order to create a pseudo- 
experimental spectrum, we have convoluted our data 
with a typical energy resolution function with a linewidth 
of ~ 1.4 meV, as measured on the x-ray inelastic scatter¬ 
ing beamline ID28 of the European Synchrotron Radia¬ 
tion Source (ESRF). The convolution gives rise to a nice 
single peak for low momentum transfers, while at high q 
we clearly see an envelope of many modes (see Fig. ITUl) . 
These pseudo-experimental data have then been fitted 
using a sum of Lorentzian functions as a model for the 
inelastic excitation. The model is then convoluted with 
the energy resolution prior to be fitted to the data. We 
find that the most intense peak of the envelope at high 
q is indeed the same mode as the single peak at low q, 
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FIG. 8: Dispersion curve for the amorphous sample ((a): longitudinal, (b): transverse) and CinASO sample ((c): longitudinal, 
(d): transverse) plotted from the dynamical structure factor. The redness of a point indicates its intensity in the dynamical 
structure factor normalized such that the most intense point for each fixed frequency is in red. 
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FIG. 9: Acoustic dispersion obtained as the maximum of the 
dynamical structure factor for different q (crosses) and from 
the ’’pseudo-experimental” method (circles) for the Amor¬ 
phous sample (longitudinal in red, transverse in purple) and 
CinASO sample (longitudinal in blue, transverse in cyan). 


giving rise to the acoustic dispersion. The acoustic dis¬ 
persions for longitudinal and transverse modes obtained 
by the two methods are in very good agreement (Fig. |9|). 
In Fig. [TT] we report the full width at half maximum 
(FWHM) of the most intense Lorentzian peak for the lon¬ 
gitudinal and transverse modes for the amorphous and 
the partially crystalline samples. We can notice that the 


linewidths obtained through the ’’pseudo-experimental” 
procedure are approximately 10 times smaller than the 
ones obtained through the mean deviation formula. This 
is expected since in the former case the fit only takes into 
account a single peak and its width. On the other hand 
the standard deviation looks at the distribution for all 
frequencies, including high frequencies where significant 
values of the dynamical structure factor can be seen. The 
linewidth as obtained from the pseudo-experimental data 
is thus more meaningfull. Beside this, the two methods 
show the same evolution for the linewidths with higher 
values for the transverse than for the longitudinal modes 
above g = 0.7 A . Once the dynamical structure factor 
is fitted, the width of the peak Aw (called damping pa¬ 
rameter in the DHO model) can be related to the lifetime 
(r*) of the phonon at frequency w. If we set r*(w) ~ 
a mean free path (MFP) can be computed as 

r(w)=Ug(w)r* (9) 

where Ug(w) is the group velocity at frequency w as far 
as this velocity is well defined (Fig. [T2|). In amorphous 
materials, eigenmodes have limited lifetime due to differ¬ 
ent contributions like scattering by the disordered orga¬ 
nization of the particles. It should be underlined that 
the dynamical structure factor cannot be fitted correctly 
for all frequencies: the loffe-Regel criterion, defined as 
hn = A/2 = Tr/q usually sets the limit between a prcm- 
agative and diffusive transport of energy and heat [^. 
This sets a critical q beyond which the phonon’s wavevec- 
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FIG. 10: Result of the convolution of the calculated S{q,ijj) 
with an experimental energy resolution function (see text for 
details) at two different wavevectors: q = 2.09 A (in red) 
and g = 4.18 A (in blue). Notice that at high q the convo¬ 
lution clearly shows an envelope of several modes when at low 
wavevector the dynamical structure factor is almost unimodal. 
The dashed lines show the peaks used to fit the spectrum. 


tor is ill-defined and the phonon cannot be seen as a plane 
wave. Moreover, Vg is ill-defined for q > 0.7 A , where 
the modes become non dispersive, thus the MFP cannot 
be calculated with Eq.[9| but other methods based on dif- 
fusivity are required j3ll |. Interestingly, it is shown from 
the dispersion relation (Fig. [5]), that the frequencies cor¬ 
responding to this wavevector correspond as well to the 
frequencies that bound the Boson peak seen in Fig. [21 
2.5 THz being the corresponding frequency of transverse 
waves and 6 THz being the frequency of longitudinal 
waves. In this frequency range, the longitudinal MFP, 
which is still well defined, decays monotonously (Fig. (121) . 
in agreement with similar calculations realized at lower 
frequencies by Larkin et al. in Ref. 31|. This decay is 
expected since the waves should have MFPs of the or¬ 
der of the size of the system at high wavelengths and be 
more and more scattered as the wavelengths tend toward 
the size of the local disorder where the wave is strongly 
scattered. 


FIG. 11: Linewidth given by the standard deviation of the dy¬ 
namic structure factor (formula [8| for the amorphous sample 
(longitudinale in red, transverse in purple) and the CinASO 
sample (longitudinale in blue, transverse in cyan). Inset: 
linewidths given by a fit of the dynamic structure factor with 
the lorentzian function. 
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FIG. 12: Longitudinal modes MFPs for the amorphous 
(in red) and CinASO (in blue) samples obtained from the 
linewidth of the dynamical structure factor with the two meth¬ 
ods (see text): pseudo-experimental (full circles) and standard 
deviation (empty circles). 


Propagation of a Wave Packet 

An alternative way to measure MFPs is to look at the 
real space propagation of a wave packet. To do this, we 
impose a wave packet to a central layer (4 A wide) of 
a sample during a MD run as done in reference Q: the 
wave packet is created by imposing a displacement to the 
particles: 


u(t) = Aexp(—^) sin(a;t)e 
2r^ 


( 10 ) 


with A = 49.10“^/a; in order to limit the kinetic en¬ 
ergy, with a maximum displacement of 49.10“^ A for 
oj = 1 THz, ensuring that anharmonic effects are negli¬ 
gible. The parameter r controls the decay of the ampli¬ 
tude of the wave with time and was set to 1 ps. This value 
is one tenth of the total simulation time of 10 ps. Lon¬ 
gitudinal and transverse wave packets can be considered 
separately depending on the relative orientation of e and 
Sx (e = Cx for longitudinal, and e.Cx = 0 for transverse 
waves). 

From the propagation of the displacement and thus 
energy along the x-axis of the sample, MFPs l(w) as a 
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FIG. 13: Linear-log representation of the kinetic energy envelopes resulting from the propagation of longitudinal waves in the 
amorphous (on the left) and ClnASO (on the right) sample. The filled symbols (limited to 20 A) are the results obtained on 
the small systems described by the dynamical matrix and the empty symbols are the results obtained on the long systems 
described by the SW potential. On the figure only some of the frequencies have been reported: 1 THz (red square), 5THz 
(purple triangle), 7.5 THz (cyan diamond), 10 THz (brown square), 12.5 THz (red circle) and 15 THz (green delta). On the 
left figure, the red line represents an exponential fit from which the MFPs are obtained, here at frequencies 1 THz. The brown 
line is a power law dehned as ajx^'^ that hts the propagation of a wave of frequency 10 THz. The purple line on the right 
figure is a guide for the eye. 


function of the frequency were obtained. For that, we 
computed first the total kinetic energy evolution with 
respect to the distance from the wave imposition layer 
{x = 0): 

N 

E{x,t) = Y^E^{t) S{\xi-x\) (11) 

i 

where Ef it the kinetic energy of the ith particle. Plot¬ 
ting E(x, t) for all the time steps during the propagation 
of the wave, we obtain a figure representing the repar¬ 
tition of the kinetic energy with time along the sample, 
Fig. [15] The initial systems are 60 A large with periodic 
boundary conditions. Once a wave reaches the border of 
a system (after 30 A), it comes back the other way and its 
energy superimposes with the energy of the leaving part 
of the wave. The simulations lasting for 10 ps, the lon¬ 
gitudinal wave has time to travel 77 A during the time 
of the simulation if we refer to the sound speed found 
previously. This distance is bigger than the size of the 
system. It limits the good fitting of the energy progres¬ 
sion in the systems to approximately 30 A, especially at 
low frequencies where the MFPs are expected to be of 
the order of the system size. To overcome this limitation 
we looked at the propagation of waves in longer systems 
made out of the repetition in the x direction of our origi¬ 
nal sample. In these longer systems, we chose to describe 
the interaction between particles using the SW potential 
that allows to take into account anharmonic effects. The 
resulting kinetic energy envelopes for both the small sys¬ 
tems described by the dynamical matrix and the big sys¬ 
tems described by the SW potential are shown in Fig.lTSl 
confirming that anharmonic effects are here negligible. 


While the energy repartition in the amorphous sample 
shows a continuous behavior, an original feature is ob¬ 
served in Fig. [T3|for CinA30. Indeed, from 15 to 50 A, 
the energy envelopes show a split in the repartition of en¬ 
ergy and this at all frequencies. In order to understand 
this split, we looked at the shape of the kinetic energy 
profile against time during the propagation of the wave, 
shown in Fig. [TH We see here that, during the propa¬ 
gation, a shoulder appears in the spatial dependence of 
the kinetic energy, giving rise to two local maxima in the 
spatial distribution of kinetic energy and consequently to 
a split in the energy envelope. To understand the origin 
of this phenomenon, it is important to look at the ge¬ 
ometry of the system. The CinA30 sample is 60 A long 
and contains a spherical inclusion in its center of radius 
24.92 A. As said, this configuration has been repeated 
several times in the x direction. Thus, in the direction of 
propagation of the wave, interfaces between crystal and 
amorphous material exist at x=5.08 A, 54.92 A and for 
the following cell at 65.08 A and 114.92 A. The interfaces 
are not plane surfaces since the inclusions are spherical 
and thus these positions are the positions of the begin¬ 
ning and the end of the crystalline inclusions. The split¬ 
ting we observe occurs exactly at the coexistence of the 
asperity and of the random matrix in the direction of 
propagation. It means that the presence of the interface 
induces a difference of speed between the high and low 
energy propagation, the former being slowed down in a 
stronger manner. This is the signature of dispersion, that 
is the existence of different velocities for different energy 
components of the wave. This could be schematically un¬ 
derstood as a plane longitudinal wave that breaks up at 
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the interface into several components with different veloc¬ 
ities similarly to what has been reported in ref. 32| in the 
case of a traveling wave packet at the interface between 
crystalline grains. Moreover, the inclusion acts as a low 
pass filter forcing the high energies (more than 1/lOth 
of the maximum kinetic energy) to pass around and not 
through it. It means that not only the frequency, but 
also the initial conditions (amplitude of the wave packet) 
affect the propagation of the wave in the composite sys¬ 
tem. It is worth noticing that the induced heterogeneity 
of the propagation is not reflected in the average MFP. 

Let’s now look at the resulting apparent MFP. At 
low frequencies, both samples show an exponential de¬ 
cay of the kinetic energy (straight red fit in Fig. [T3l) . At 
higher frequencies and large distances, we can observe a 
transition to a power law behavior, for the longitudinal 
wave packet only (e = Cx)- The characteristic frequency 
that marks the separation between these two behaviors 
is 2.5 THz. This change of trend can be interpreted as 
the crossing from simple to multiple scattering [33| . How¬ 
ever, multiple scattering behavior is usually fitted by an 
Ohm’s like law {l/x) [33| and we had to use a power law 
with a Ri 3 for w > 2.5THz see Table HIl) to match 
the energy repartition of our data. Concerning the low 
frequency behavior of the transverse wave packets (with 
e.Bx =0), it can be fitted with a decreasing exponential 
function of the form B x exp {—x/l) known as the Beer- 
Lambert law, from which I, the MFP of the wave, can be 
found. For the CinASO, since the energy envelope is not 
homogeneous, we chose to fit the average of the curve 
in order to obtain the MFPs. These l{uj) can be com¬ 
pared to those found through the fitting of the dynamical 
structure factor, for q < 0.7 A . The behavior is 
similar, although the mean-free path obtained from the 
standard deviation formula of the structure factor is un¬ 
derestimated. For the longitudinal waves, an almost con¬ 
stant MFP can be observed in all cases at low frequency 
before a decrease starting at 2.5 THz with a minimum 
at 7.5 THz, followed by a regain around 9 THz already 
observed in Ref. 3^. The longitudinal MFP obtained 
from the wave-packet propagation decays from 20 A to 
2 A. The low frequency value is probably limited by the 
finite box size: indeed it is much smaller than the one 
obtained from the linewidth which is in agreement with 
Ref. 3^. The MFP for the transverse waves is always 
smaller than the one for the longitudinal waves, start¬ 
ing at 6 A at small frequencies down to 0.6 A at larger 
frequencies. It follows a similar trend with a saturation 
followed by a regain at 9 THz. This regain is unexpected 
since the transverse waves have been shown to disappear 
after the boson peak [2^ but could be explained only by a 
coupling between longitudinal and transverse modes dur¬ 
ing the dynamical propagation of the wavepacket. The 
largest discrepancy between amorphous and CinASO sam¬ 
ple occurs for transverse waves in the 2 to 12 THz range. 
From these sets of data however it is not possible to defini- 


Lo (THz) 

0.3 to 1 

1.5 

2 

2.5 

3 

3.5 

5 

q: 

- 

1.9 

2 

2.2 

2.5 

2.8 

3 

OJ (THz) 

6 

7 

9 

10 

12 

15 

20 

a 

3 

3 

3.1 

3.2 

3.3 

3.2 

2.9 


TABLE II: Exponent a of the power-law fit of the large dis¬ 
tance energy envelope, as a function of the frequency lj, for 
longitudinal wave packets. 



x{A) 

FIG. 14: Kinetic energy repartition for different time steps 
(t=n*800 fs) as a wave of frequency 5 THz propagates in the 
CinASO sample. 


tively assess which sample between the amorphous and 
CinASO is more resistive to the propagation of waves and 
energy, because the effect of the inclusion on the average 
mean-free path appears to be very small. 

While MFP obtained from the wave propagation and 
from the Lorentzian fits are similar, the detailed study 
of the waves propagation show evidence of spatial het¬ 
erogeneities. In the amorphous sample, we have ob¬ 
served a transition from simple scattering-like to multi¬ 
ple scattering-like in the attenuation of energy of a lon¬ 
gitudinal wave-packet, at a scale comparable to the mea¬ 
sured MFP. In the presence of the inclusion (whose size 
is larger than the MFP of the amorphous matrix), the 
crystallite acts as a low pass filter for the kinetic energy, 
adding fluctuations to the spatial envelope of the energy 
at larger scale than the MFP. These heterogeneities can¬ 
not be catched by the MFP, but still they must play 
an important role in thermal transport, thus raising the 
question of the correct modeling of thermal transport in 
composite systems. 
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FIG. 15: MFP of transverse (squares) and longitudinal 
(crosses) plane waves as a function of the frequency for the 
amorphous (in red) and CinA30 (in blue) samples obtained 
through the propagation of wave packets. 


CONCLUSIONS AND DISCUSSION 


In this paper, we studied in detail the role of nano¬ 
metric crystalline inclusions on the vibrational properties 
of amorphous samples (a numerical model of Si sample). 
We find that the VDOS of a sample with inclusions is 
approximately described by the sum of the VDOS of the 
amorphous and of the crystalline samples, although sub¬ 
tle differences still persist. However, the coefficients of 
the linear combination are not proportional to the volu- 
mic fraction of crystal in the sample, indicating that this 
is not the right parameter. The main effect of the inclu¬ 
sions on the density of states is to lower the number of 
resonant modes at frequencies close to 12.5 THz, a fre¬ 
quency where the density of vibrational modes is very 
low in crystalline samples. The average MFP measured 
here in different ways appears to be insensitive to this 
peculiar frequency, and does not show a strong depen¬ 
dence on the inclusion. One visible effect of the presence 
of the inclusion is a decrease of the MFP of transverse 
waves for frequencies between 2.5 and 10 THz. These fre¬ 
quencies are close to the frequencies bounding the Boson 
Peak in the amorphous sample, that were related in this 
paper to the frequencies at which transverse (2.5 THz) 
and longitudinal (7.5 THz) modes are significantly scat¬ 
tered by the amorphous structure. Interestingly, from a 
close inspection of the MFP as obtained looking to the 
propagation of wave packets, we find that longitudinal 
wavepackets are already multiply scattered at large dis¬ 
tance for frequencies smaller than 7.5 THz. However the 
related longitudinal mean-free path is far larger than the 
transverse one, and no effect is visible in the dynamical 
structure factor. The two frequencies mentioned above 
appear in the respective dispersion relations for the same 


value of the wavevector q = 0.7 A ^ that is a wavelength 
\ = 2tt/ q = 0 K. This observation rises the question of 
the existence of a characteristic length for waves scatter¬ 
ing in amorphous samples 23|. This mesoscopic length 
is larger than the typical length scale for plastic rear- 
angement (5 to 6 A [1^) but of the same order of the 
mean-free paths and of the distance above which the Beer- 
Lambert law for acoustic attenuation breaks down to a 
Ohm’s like law. The fact that the transition from simple 
to multiple scattering is seen only in longitudinal waves 
could be due to the fact that the transverse mean-free 
path is always smaller that this characteristic distance. 
On the other hand, this length cannot be linked easily to 
any structural feature in our amorphous samples: larger 
than the average first-neighbour distance, it could be re¬ 
lated to long range correlations due to disorder and 
elastic heterogeneities [s^. However, further work is still 
needed to understand if it can be related to some charac¬ 
teristic local order existing in amorphous silicon. 

The main effect of the crystalline inclusion is to act 
as a low-pass filter for the kinetic energy. At the inter¬ 
face between amorphous and crystalline parts, the high 
kinetic energies are slowed down more strongly than the 
low kinetic energies. This results in an amplified hetero¬ 
geneous waves propagation, and in enhanced fluctuations 
in the energy envelope, thus raising the question of the 
mean-free path definition in this case, as well as its role 
on thermal transport 36|. A deeper understanding of 


the microscopic origin of this unexpected kinetic effect 
would allow to design nanostructured samples with very 
high acoustic attenuation. 

More work is needed for investigating the effect on 
thermal transport of size and shape of the crystal 
inclusion as well as the possible interaction when there 
are more than one inclusion. Due to the relativily small 
size of the systems that can be studied using MD, other 
methods are advisable as it has been done by Zhang et 
al. in Ref. 371. 
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